Effective Dynamics in Bianchi Type II Loop Quantum Cosmology 

Alejandro Corichi 1 ' 2 '^ and Edison Montoya 3 ' 

1 Centro de Ciencias Matemdticas, Universidad National Autonoma de Mexico, 
UNAM-Campus Morelia, A. Postal 61-3, Morelia, Michoacdn 58090, Mexico 
2 Center for Fundamental Theory, Institute for Gravitation and the Cosmos, 
Pennsylvania State University, University Park PA 16802, USA 
3 Instituto de Fisica y Matemdticas, Universidad Michoacana 
de San Nicolas de Hidalgo, Morelia, Michoacdn, Mexico 

We numerically investigate the solutions to the effective equations of the Bianchi 
II model within the "improved" Loop Quantum Cosmology (LQC) dynamics. The 
matter source is a massless scalar field. We perform a systematic study of the space 
of solutions, and focus on the behavior of several geometrical observables. We show 
that the big-bang singularity is replaced by a bounce and the point-like singularities 
do not saturate the energy density bound. There are up to three directional bounces 
in the scale factors, one global bounce in the expansion, the shear presents up to four 
local maxima and can be zero at the bounce. This allows for solutions with density 
larger than the maximal density for the isotropic and Bianchi I cases. The asymp- 
totic behavior is shown to behave like that of a Bianchi I model, and the effective 
solutions connect anisotropic solutions even when the shear is zero at the bounce. 
All known facts of Bianchi I are reproduced. In the "vacuum limit" , solutions are 
such that almost all the dynamics is due to the anisotropies. Since Bianchi II plays 
an important role in the Bianchi IX model and the Belinskii, Khalatnikov, Lifshitz 
(BKL) conjecture, our results can provide an intuitive understanding of the behavior 
in the vicinity of general space-like singularities, when loop-geometric corrections are 
present. 

PACS numbers: 04.60.Pp, 04.60.Bc 



I. INTRODUCTION 



Loop Quantum Gravity (LQG) [T] has recently emerged as a strong candidate for a 
quantum theory of gravity. One of the main motivations for such a theory is to provide 
a solution to the big-bang singularity. Even when the full theory has still little to say 
about the initial singularity, a symmetry reduced theory, namely Loop Quantum Cosmology 
(LQC) [2HI], has been extremely successful at providing precise answers to that question. 
The theory is constructed by applying the methods of LQG to a symmetry reduced sector 
of general relativity. As examples of these reduced configurations, several authors have 
studied cosmological models with a massless scalar field and geometrically isotropic [5]- 
[9], homogeneous and anisotropic [T0HT3] . and some inhomogeneous cosmologies [14] . The 
common theme among these models is that they resolve the big bang singularity. The way 
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singularity resolution occurs is by means of physical observables (in the sense of Dirac) whose 
expectation values (or spectrum) have been shown to be bounded [ToT - flS] . These results 
benefit from uniqueness results that guaranty the consistency of the so called "improved 
dynamics" [IBl US]. 

In a sense, isotropic LQC can be seen as a realization of one of the main objective of 
LQG, namely, to solve the big-bang singularity. In this case the singularity is replaced 
by a bounce that occurs precisely when the matter density enters the Planck regime. At 
this energy density, the quantum effects create a repulsive force, the would-be singularity is 
avoided and the resulting 'quantum' spacetime is larger than one might be led to believe. 
As has been shown in detail, when the density decreases, the state very quickly leaves this 
quantum regime, and the universe returns to being well described by general relativity [2]. 
With the final goal in mind of investigating the most general issue of singularity resolution 
within LQG, one expects to gain useful insights from results obtained for less symmetric 
models. The simplest such model is given by Bianchi I cosmologies [H] where no spatial 
curvature exists. Even for this case we do not possess yet a full exact evolution of the 
quantum equations of motion. 

For both isotropic and anisotropic models, a very useful tool has been the use of the 
so-called effective description, a 'classical theory' (in the sense that it does not contain h) 
that has information about the geometric discreetness contained in loop quantum gravity. 
In isotropic models, the solutions to the effective equations have been shown to approximate 
the dynamics of semiclassical states in the full quantum theory with very good accuracy [7]- 
IH 120] • In the case of anisotropic models we also expect that the effective solutions play an 
important role for describing the evolution of semiclassical states. Thus, we shall adopt the 
viewpoint that it is justified to study the effective dynamics for those cosmological models, 
as a way to learn about the full 'loopy' quantum dynamics. 

So far, the only anisotropic model that has been explored in detail, through the so called 
effective equations, is the LQC Bianchi I model [ID] . It is then natural to explore this issue 
for more complicated anisotropic models. In this paper we will study the effective equations 
obtained from the "improved" LQC dynamics of the Bianchi II model |12| . The Bianchi II 
cosmological model represents the simplest case that posses spatial curvature, from which 
the Bianchi I model can be recovered when a parameter measuring the spatial curvature 
contribution is 'switched off'. The Bianchi II model possesses another interesting feature, 
namely, it lies at the heart of the Belinskii, Khalatnikov, Lifshitz (BKL) conjecture [2TH23] . 
which suggest that, as one approaches space-like singularities, the behavior of the system 
undergoes Bianchi I phases with Bianchi II transitions. One question that remains open 
though is whether this BKL behavior will survive in the effective theory. That is, will the 
oscillations between Bianchi I phases occur far from the Planck scale? or will the loop- 
geometric effects prevent this 'mixmaster' behavior to manifest itself? From this point of 
view, it is important to study in detail the Bianchi II solutions. The most natural strategy to 
gain this intuition is to perform a systematic study of the solutions to the effective equations, 
under the assumption that they describe correctly the quantum dynamics. But even if the 
effective solutions do not describe correctly the quantum dynamics of the semiclassical states 
in some regime, we need to study them in detail in order to compare them to the full quantum 
dynamics -when available- and prove their validity. One can also expect that this study 
will shed some light on the larger issue of understanding generic space-like singularities in 
LQG. 

The purpose of this paper is to study in a systematic way the space of solutions of the 
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effective equations. The matter source that we shall consider is a massless scalar field that 
plays the role of internal time. The objective is to understand the singularity resolution, 
the asymptotic behavior and the relation between the Bianchi II and Bianchi I models. The 
strategy will be to take limiting cases and compare them with known solutions. This will 
allow us to understand the new insights of the Bianchi II model. All the information will 
come from the set of observables that we define, namely, directional scale factors, Hubble 
parameters, expansion, matter density, density parameter, shear, shear parameter, Ricci 
scalar, curvature parameter and Kasner exponents. With these tools in hand we shall 
compare the classical and the effective solutions. Later on, the isotropic limit will offer 
interesting new insights into the Bianchi II dynamics, while the Bianchi I limit will allow us 
to confirm that our solutions agree with the previous results of [10]. The symmetry reduction 
to the Locally Rotationally Symmetric (LRS) Bianchi II model will give us a way to find 
generic solutions with maximal density, at the bounce, larger than the critical density derived 
in the isotropic case [7J [T7j. Further, the vacuum limit will allows us to study the extreme 
solutions where all the dynamical contributions come from the anisotropies. Finally, due 
to the fact that our work is numerical we will show that our solutions converge and evolve 
on the constraint surface. The convergence is an important issue that needs to be shown, 
because one needs to ensure that the numerical methods are well implemented, and that the 
numerical solutions converge to the analytical solutions. 

The structure of the paper is as follows. In Sec. II we recall the classical theory in metric 
and connection variables, together with the basic observables to be studied. In Sec. Ill we 
introduce the effective theory and compute its equations of motion. Numerical solutions are 
explored in Sec. IV, where they are systematically studied taking the classical, isotropic, 
Bianchi I, maximal density and vacuum limits. We end with a discussion in Sec. V. There 
is one Appendix where we study the convergence and evolution of the conserved quantities. 
Throughout the manuscript we assume units where c = 1, the other constants are written 
explicitly. 



II. CLASSICAL DYNAMICS 

The metric for a Bianchi II model can be written as 

ds 2 = -A(r) 2 dr 2 + ai (r) 2 (da; - az dyf + a 2 (r) 2 dy 2 + a 3 (r) 2 dz 2 , (2.1) 

where the parameter a allows us to distinguish between Bianchi I (a = 0) and Bianchi 
II (a = 1) cases. Bianchi I cosmological solutions are also interesting, since they give 
information about the asymptotic behavior of Bianchi II. Classically, the Bianchi I case 
with a massless scalar field has solutions given by [23], 

ds 2 = _ df 2 + t 2k, dx 2 + t 2k 2 dy 2 + f 2k 3 ^2 ^ 

the so-called Jacobs stiff perfect fluid solutions. 1 Here, the parameters ki are known as 
the Kasner exponents satisfying k\ + k 2 + k 3 = ±1 (the minus sign is inserted, for future 
reference, as we want to take into account the change of direction of the expansion at the 



1 The metric form (2.2) is taken from [25], where one can find a large class of analytical solutions to spatially 
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bounce) and kf + k 2 + A;f + k\ = 1. In the literature, a massless scalar field is also known as 
"stiff matter" , 2 and satisfies the equation of state P = p, where P is pressure and p is the 
energy density. All these solutions have an initial singularity at t — 0, that can be of four 
types [21]: 

1. Point type singularity. This means that 01,02,03 — > as t — > 0. This happen in the 
Jacobs solutions when k\, k 2 , &3 > 0. 

2. Cigar type singularity. These occur when 01,02 — » and 03 — > 00 as t — > 0. This 
happens when fci, k 2 > and fc 3 < 0. (cyclic on 1,2,3) 

3. Barrel type singularity. Defined by 01,02 —> and 03 approaching a finite value as 
t — > 0. This happen when A^, fc 2 > and A>3 = 0. (cyclic on 1,2,3) 

4. Pancake type singularity. In this case, ai — > and a 2 ,a 3 approaches a finite value as 
£ — > (cyclic on 1,2,3). This kind of singularity is not realized in the Jacobs solutions 
(nor the Bianchi II with a massless scalar field, or "stiff" matter) (26]. This can be 
understood easily. This singularity happens when k\ > and k 2 = k 3 = (cyclic on 
1,2,3) in the Bianchi I case, satisfying the constraint equations: ki + k 2 + k% = 1 that 
gives k\ = 1. Now, the condition kf + k\ + fcf + k 2 ^ = 1 implies that k<p = 0, which 
means that there is no matter. Thus, this type of singularity is not possible if there is 
stiff matter present. 

These names for the different types of singularities, introduced by Thorne in |27J, refer to 
the change of the shape of a spherical element as the singularity is approached. The Jacobs 
solutions will be useful in our analysis, since Bianchi II is past and future asymptotic to 
Bianchi I 3 (see, for instance [281 12H] and section 9.3 of [30]). Additionally, Bianchi II is a 
limiting case for the effective equations of Bianchi II that come from its quantization [12] 
in the loop quantum cosmology (LQC) framework. Then, in the classical region, Bianchi I 
(Jacobs) solutions are limiting cases for the effective Bianchi II that arises from LQC. 

We will now rewrite the classical theory in terms of triads and connections in order to 
connect with the effective theory that comes from the quantum theory. To do this we use 
the fiducial triads and co-triads and introduce a convenient parametrization of the phase 
space variables, Ef, A % a given by 

E? = p^V^^m and A\ = c t L^ ^ (2.3) 

without sum over i, where Vo = LiL 2 L 3 is the fiducial volume and Lj the fiducial lengths 
with respect to the fiducial metric °g ab := S^ld^l with co-triads 

°u;i = (da;) a -az(d 2 /) a , °co 2 a = (dy) a °ul = (dz) a (2.4) 

and triads 

*-(!)'■ (=)" + (£)"- °^©"< (2 - 5) 

2 In the original article [53], this solution is called "Zel'dovich universe". 

3 This is true when the matter is a massless scalar field. For other matter the asymptotic solutions are 
different. 
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with Lie bracket [ e2, e3] = — °ei, 4 , with a = 1. A point in the phase space is now coor- 
dinatized by eight real numbers (pj, q, 0, p^), with (/> the scalar field and p^ its conjugate 
momentum. The Poisson brackets are given by 



{ci, pj} = 8nG-f6i 



{0,M 



(2.6) 



where 7 is the Barbero-Immirzi parameter. The Hamiltonian formulation will be complete 
with the Hamiltonian constraint, which for a lapse function N = \/\p\P2P3\ reads 



C 



H 



8nGr 



P1P2C1C2 +P2P3C2C3 +PiP3Cic 3 +aep 2 p3Ci - 



(2.7) 



where again a distinguish between Bianchi I and Bianchi II, e = ±1 depending on whether 
the frame e\ is right or left handed (in our solutions we assume e = 1, i.e. pi > 0). The 
equations of motion are given by the Poisson brackets with the Hamiltonian constraint 

(2.8) 
(2.9) 
(2.10) 

(2.11) 



(2.12) 

(2.13) 

(2.14) 
(2.15) 



Pi 


= T 1 {piP2C2 + P1P3C3 + aep 2 p 3 ), 






P-2 


= l~ 1 (P2PiCi +P2P3C3), 






P3 


= l~ 1 (P3PlCl+P3P2C2), 








= -7- 1 


P2C1C2 + P3C1C3 + (l + 7 2 )( 
L 2p 1 


aP2P3\2~ 

Pl 


) 




= -7- 1 


r 1 

P1C2C1 + P3C2C3 + otepzCi - — 
L 2p 2 


p 1 J 


C3 


= -7" 1 


r 1 
P1C3C1 + P2C3C2 + aep 2 c 1 - — 
L 2p 3 


p 1 J 


Pcf, 


= 


=^ P<i> = constant, 









= ~P<t> 


<P = ~P<P T > 







where r is called the harmonic time (with lapse iV = W \piP2Ps\)- The last equation shows 
that the field can be used as internal time. The harmonic time r is related to the cosmic 
time t (with lapse N = 1) by the equation 



d 

d7 



(2.16) 



^/\PlP2P3\ dr 

It is with respect to this time that we will define the observable quantities. The derivative 



respect to the cosmic time will be denoted by dO/dt = O', then O' = 0/ y \piP2P3\- From 
the equations of motion it is straightforward to show that the classical solutions posses the 
constants of motion 



cm + C2P2 

cm + c 3 p 3 

C3P3 - C2P2 



"12 
"13 
"32 



"13 - «12 



(2.17) 
(2.18) 
(2.19) 



4 We are using the notation from [T^] and chapter 11 of [31]. Note that this is not the typical choice 
([°e2,°e3] = °e%) for Bianchi II. This choice only implies a change from one invariant set (rii = 1) of 
Bianchi II to the other one (ni = —1), but the physical properties are the same, given that the equations 
of motion have this discrete symmetry. For more details see [30H32) . 
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with ai2,«i3 constants. This will allow us to solve analytically equations (2.9) and (2.10), 
t>2 = l~ 1 P2(PiCi +P3C3) = 7 _1 P2ai3 P2= P° ex P y^J ' ( 2 - 20 ) 

P3 = 7"V3(PlCl + P2C2) = 7"V3ttl2 P3=P°3 exp " ( 2 ' 21 ) 

The existence of this exact solutions gives us the opportunity to compare the numerical 
solutions with the analytical ones for p 2 and p%. Also, we can check that during the 
evolution, for the numerical solutions, a,\ 2 and (X13 remain constant. 

In order to determine how the classical singularities are resolved and how the effective 
equations evolve, the quantities that we will study are: 

1. The directional scale factors 



L -i IPjPk 



Pt 

with i ^ j / k ^ i and Pi > 0, i — 1, 2, 3. This makes easy the comparison between 
the classical and the effective solutions. 

2. The directional Hubble parameters 

Hi= < = 1 (Pj L + P'k_P'i 



CLi 2 \Pj p k Pi 

with i 7^ j / k 7^ i. This quantity tells us when each direction bounces (Hi = 0) or if 
these directions are contracting (Hi < 0) or expanding (Hi > 0). 

3. The expansion 

V 

8 = — = Hi + H 2 + H 3 . 

This quantity gives the total expansion rate and determines when there is a global 
bounce (9 = 0). An equivalent quantity is the mean Hubble parameter 

H= 6 - = ^(H 1 + H 2 + H 3 ). 

4. The matter density 

P\ P\ 

P ~ 



2V 2 2p lP2 p 3 - 

If the singularities are resolved then this quantity must be finite thoughout the evolu- 
tion. Other quantity that measures the dynamical importance of the matter content 
is the density parameter Q, defined by 

This parameter is related to the Kasner exponents in Bianchi I by the equation Q = 
\k\ = |(1 - kf - k\ — kf) (see, for instance [2"4]). 
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5. The shear 

1 3 1 

a 2 = a ab a ab = i[(tf a - H 2 ) 2 + (H, - H 3 ) 2 + (H 2 - H 3 ) 2 ] = ^H 2 - -9 2 . 

i=i 

Note that this definition of cr 2 differs from the standard definition a 2 = \cr a b(J ab . 
Another important quantity is the shear parameter 

2 _3a 2 a 2 

that measures the rate of shear (i.e. anisotropy) in terms of the expansion. In Bianchi 
I the expansion satisfies 9 = t -1 (using the fact that V = a\a 2 a 3 = t kl+k2+ki = t, 
without putting explicitly the units). Then, the shear parameter reduces to the relation 
£ 2 = \t 2 o 2 = \V 2 a 2 = \d°o 2 = 9S|j, where E| r := \a 2 a 6 was the shear parameter 
used in Bianchi I [101 [TTJ, with a := (a^a^ 1 / 3 the mean scale factor. 



6. The Ricci scalar for the Bianchi II metric, Eq. (2.1), is given by 



CL-^ ^3 CL-^CL<2 (X-^CL^ ^2 ^3 \ 2 ^ ^ 



^^H^^^^^ -oc 2 -^. (2.24) 
\ai a 2 a 3 a\a 2 a\a 3 a 2 a 3 J 2a 2 a 3 

When a = it reduces to the Ricci scalar for Bianchi I. In terms of the new variables 
(ci,pi) it has a simple expression 

R= p1 + pA + pA-\ x \ (2.25) 

Pi P2 P3 2 



where x = a. This equation can be rewritten in terms of other observables, 



Pi 

R = 26' + a 2 + \e 2 - \x 2 . (2.26) 

This relation provides the easiest way to calculate the Ricci scalar numerically. On 
the classical solutions 6' is given by the Raychaudhuri equation 

0' = --6 2 -a 2 -16irGp. (2.27) 

One important feature of the Bianchi II model is that the spatial curvature is different 
from zero, thu, we can introduce other quantity that give us information about the 
dynamical contribution due to the extrinsic curvature, namely the curvature parameter 
K, given by 

Our choice of Q, S 2 and K is motivated by the fact that, on the classical solutions for 
Bianchi II, they satisfy the equation 

VL + Y? + K = 1. (2.29) 



7 



These parameters have the 'problem' that they are infinity at the bounce (6 = 0) 
by definition, so they are not very useful to explore that regime. But, since we are 
interested in its asymptotic behavior for large volume, and the information they can 
give us there, this pathological behavior at the bounce is not relevant and does not 
reflect any problem with the singularity resolution. 

7. The Kasner exponents 

h = ^ (2.30) 

These parameters are very useful to determine when the solutions have a Bianchi I 
behavior. We have taken the absolute value in 9 because we want that different signs 
in Hi specify i f the directions are expanding {ki > 0) or contracting (ki < 0). In order 



to prove Eq. (2.30) we use that, in Bianchi I, = t ki Li (using explicitly the fiducial 
lenghts), then 



_ a'/- fe * _ a'it _ a\t _ a%t _ a\V _ a-V _ Hi 
Li t k *Li a,i aiV aiV aj/' 6 

with V = a\a 2 a 3 = Vot. 



(2.31) 



One important remark is that all the previous expressions to calculate the observable 
quantities apply also to the effective solutions. The only difference is the calculation of 9' in 
the effective theory, for which it is necessary to use the effective equations of motion. That 
will be shown in the next section. 

To complete the classical picture we give the relations between phase space variables and 
metric variables, which are given by 

Pi = a 2 a 3 L 2 L 3 , p 2 = aia 3 LiL 3 , p 3 = a x a 2 LiL 2 , (2.32) 

and 

a a?L? 

a = 7 £iai#i + - — VV > ( 2 - 33 ) 
2 a 2 a 3 L 2 L 3 

c 2 = ^L 2 a 2 H 2 --^ } (2.34) 
l a 3 L 3 

c 3 = ^L 3 a 3 H 3 - ^^L_i ; (2.35) 
l a 2 L 2 



assuming > 0,pi > and e = 1. Relations (2.32) are satisfied at the kinematical level 
whereas Eqs. ( 2.33|2.34|2.35 ) are satisfied at the dynamical level, i.e., on the space of 
solutions. These relations can be shown using the Hubble parameters (with explicit fiducial 
lenghts), 

H l = ^ = 1 (k + Ph_Pj) (2 .3 6 ) 
di dt 2 y /p{p 2 ~p 3 ~ \pj p k pi J 

dai _ 1 fpj + pk_ Pi\ ^ 2 3 ^ 



dt 2piLi \pj p k pi 



S 



with i ^ j ^ k ^ i. If we put the equations of motion (2.8, 2.9 2.10) into this expression, 
we get 



dai a 1 ap 2 p3 

-TT = -^ + fsign(« — — 

at ^Li 2p i L i 'j pi 



(2.38) 



where fsign(i) = — 1 if i = 1 and fsign(z) = 1 if % = 2, 3. Using the relations between ai and 
Pi, Eq. (2.32), and the definition of the Hubble parameters we get 



c t = lL~ - fsign(z)^-aLiai 
at 2pj 

t tt f - ( -\ a L i a i 

= -yLidiHi - tsign(z) 



2 LjLkCijCik 



(2.39) 
(2.40) 



These are precisely Eqs. (2.33, 2.34, 2.35). 



III. EFFECTIVE DYNAMICS 



The effective theory was derived from the loop quantization of Bianchi II defined in 
reference [12], using the procedure outlined in [33J. Taking a right-hand frame e" (i.e. 
e = 1) and the lapse function iV = y/p^pipi, the effective Hamiltonian constraint is given 
by, 



C 



H 



PlP2P3 

8ttG 7 2 A 2 
1 

+ 



sin niC\ sin p 2 C2 + sin p 2 C2 sin /i 3 C3 + sin /Z3C3 sin \i\C\ 



8ttG 7 2 



a{p2P-i) 



3/2 



sm/i lCl - (1 + 7 ) — 

Ay/Pi V 2 Pl / 



P% 



with 



A 



Pi 



fJ>2 



A 



P2 



/'3 



A 



Pa 



(3.1) 



(3.2) 



P2P3 V P1P3 V P1P2 

The value of A is chosen such that A 2 = A corresponds to the minimum eigenvalue of the 
area operator in loop quantum gravity (corresponding to an edge of "spin 1/2"). With this 
choice the free parameter becomes A = 4v / 3tt 7 ^ p1 . Since sinpjCj < 1 the matter density 

satisfies 



P 



2V 2 



Pi 



2piP2P3 



< 



+ 



87rG 7 2 A 2 8ttG 7 



x (l + 7 2 )x 2 
A 4 



with x 



a 



P2P3 
Pi 



The maximum of the expression in square brackets is attained at x 

Pmatt < 1.315p crit ~ 0.54p P1 , 



(1+7 2 )A 



(3.3) 

0.83, then 
(3.4) 



with p crit = 87rG y A 2 ~ 0.41ppi the critical density found in the isotropic case (6] and ppi = 
mpi/£ P1 is the Planck density. This shows that the matter density is bounded. This bound is 
higher than the one found in Bianchi I p mat t < 0.41ppi and the isotropic case p mat t ~ 0.41ppi. 
Furthermore, the density in all the solutions in Bianchi I with shear different from zero has a 
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bounce density less than its value in the isotropic solution. Then there is an open question: 
Are there generic solutions in which the matter density is larger that its value in the isotropic 
solutions? We will show that the answer is in the affirmative, which leaves us with another 
open question: How do we find the solutions that saturate the matter density? These kind 
of solutions are shown in section IIVDI 



If we set a = into Eq. (3.1) we recover the Hamiltonian constraint for Bianchi I [TT] . 
Also, if we take the Bianchi II case, a = 1, we can recover Bianchi I as a limiting case 
when x — > 0, or equivalently p\ p 2 p 3 (in metric variables this condition is expressed as 
a.20-3 ^> cti). It is important to have in mind that the Bianchi I model is a limiting case and 
is not contained within the Bianchi II model. 

The equations of motion for the effective theory are given by Poisson brackets with the 
Hamiltonian constraint, 



Pi 

P3 



Pi 

P\ { 
Jfl2 

P\ , . _ . - N 

— (sm/iiCi + sin// 2 c 2 J cos /i 3 c 3 , 



(sin /z 2 c 2 + SU1/X3C3 + Ax) cos faCi, 



v sin/i!Ci + sin /x 3 c 3 ) cos/i 2 c 2 , 



(3.5) 
(3.6) 
(3.7) 



Ci = ~^^2~ |^ 2(sm /iiCi sm / u 2 c 2 + sm/iiCi sm /i 3 c 3 + sin/i2C2sm / u 3 C3) 

+ /iiCi cos /2iCi (sin /i 2 c 2 + sin/2 3 c 3 ) - /i 2 c 2 cos /i 2 c 2 (sin /2iCi + sin /i 3 c 3 ) 
- /2 3 c 3 cos/2 3 c 3 (sin/i 1 Ci + sin Jx 2 c 2 ) + X 2 x 2 (l + 7 2 ) 

+\x{p,\Ci cos Ji\C\ — sin Ji\C\) 



(3.8) 



c 2 



P1P3 
2 7 A 2 



2 (sin fiiCi sin // 2 c 2 + sin \i^C\ sin /i 3 c 3 + sin /i 2 c 2 sin /i 3 c 3 ) 



/iiCi cos /iiCi (sin /i 2 c 2 + sin /i 3 c 3 ) + /i 2 c 2 cos /i 2 c 2 (sin [i x c x + sin/^ 3 c 3 ) 
- / u 3 c 3 cos/i 3 c 3 (sin/iiCi + sin /i 2 c 2 )] - A 2 x 2 (l + 7 2 

\x{ji\C\ cos Ji\C\ — 3 sin Ji\C\) 



(3.9) 



P\P2 



c 3 = -7^^ 2(sin/iiCi sin /x 2 c 2 + sin /iiCi sin /x 3 c 3 + sin /i 2 c 2 sin /x 3 c 3 ) 

- /iiCi cos /iiCi (sin /i 2 c 2 + sin/2 3 c 3 ) - /i 2 c 2 cos/i 2 c 2 (sin/2 1 Ci + sin/2 3 c 3 ) 
+/2 3 c 3 cos /Z 3 c 3 (sin \x x c x + sin /2 2 c 2 )] - A 2 x 2 (l + 7 2 ) 



— Ax(/iiCi cos /iiCi — 3 sin fi±Ci) 
Finally, for the matter we have 

= 



constant, 



(3.10) 



(3.11) 
(3.12) 
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Note that the equations for the matter part are equal to the classical ones, so the field 
also plays the role of internal time. The field momentum p^ is conserved, and its value 
coincides with the classical value. From the triad equations (p\) we can see that the bounce 
in direction i occurs when cos(/2jCj) = 0, which can be satisfied at different times for each 
direction. These assertions do not imply that there is more than one bounce in the matter 
density; the density only bounces one time, which is when the expansion 9 is zero. This is 
called the global bounce. 

We can use the equations of motion to give an explicit formula for 9' (derivative of the 
expansion with respect to cosmic time), which is necessary to calculate the Ricci scalar. It 
is straightforward to show that 

1 f 3 

9' = — < ^2 psin/^Ci + cos fiiCiisfci' (ijCj + sin' p, k c k )] 

+Al > M!ftCl « + !l_Wl (3.13) 

2 \P2 P3 Pi J J 

where j ^ i ^ k ^ j and sin'(/2jCj) = cos(p,iCi)[p' i Ci + fii^), with = —fiiHi. It would be in- 
teresting if one could rewrite Eq. (3.13 ) in terms of observable quantities (Hi, 9, p, cr 2 ), which 



would represent a generalization of the Raychaudhuri equation. Also, from the equations 
for p\ and q we get that 

C3P3 - c 2 p 2 =: "32 (3.14) 



is conserved and its value is the same than the classical one, as given by Eq. (2.19). The 
conserved quantities (p^ and a 32 ) can be used to check that numerical solutions are evolving 
correctly on the constraint surface. 



IV. NUMERICAL SOLUTIONS 

In this section we show the numerical solutions for the Bianchi II model. These equations 
admit different limits than can be used to check the accuracy of the solutions and explore 
the new insights that Bianchi II offers. In order to systematically study the solutions we 
need to have in mind the following facts: 

1. In all the solutions, the constraint (or equivalently p<$) and CK32 = c^pz — c 2 p 2 are 
conserved quantities. 

2. In the classical solutions a% 2 and are conserved quantities. 

3. When a = or x = \/ PiVilvX ~^ (with a = 1), Bianchi II reduces to Bianchi I. 

4. When p\ = p 2 = p% and c\ = c 2 = C3, Bianchi I reduces to the isotropic case. 

5. In the isotropic case all the solutions to the effective equations have a maximal density 
equal to the critical density p crit = 87rG 3 y2A 2 ~ 0.41p P1 . 

6. In the Bianchi I limit we will expect a maximal density less than the critical density. 
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FIG. 1: Comparison between classical and effective solutions. Here it can be appreciated that near 
the big bang the classical solutions have an infinite density because the scale factors go to zero. 

7. Classical and effective solutions must be equal far away from the bounce. In fact 
the effective solutions to Bianchi I must connect two classical solutions with Kasner 
exponents related by 



as was shown by Choiu [ID] . 

8. In the Bianchi II model we expect a maximal density less than 1.315p cr it. 

9. The classical solutions diverge. 

10. The numerical solutions must converge. 

Using these facts we shall now explore the Bianchi II solutions. We start from the classical 
limit showing that the effective solutions have a bounce and reduce to the classical ones far 
away from the bounce. Next, we explore the isotropic limit included into the Bianchi I 
limit when there are no anisotropics. Later on, we shall add anisotropics to the Bianchi I 
limit and show that they reproduce the known solutions [10J. Then we pass to the Locally 
Rotationally Symmetric (LRS) model of Bianchi II and explore how to find the solutions 
with maximal density at the bounce. Finally we study the solutions in the vacuum limit 
with maximal shear. 

To perform the analysis, solutions are plotted as functions of the cosmic time. All the 
integrations were made using a Runge-Kutta 4 method 5 . The units are: h = l,c = 1,G = 
1,7 = 0.23753295796592, L x = 1, L 2 = 1, L 3 = 1, e = 1, the time step used is dt = 5 x KT 6 . 
In all figures the density is plotted with normalization p/p C rit- 



5 The program is available by request. 



k l7 k 2 ,k 3 -»• ki - 




(4.1) 
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A. Comparison Between Classical and Effective Solutions 



In figure [T] the density and two scale factors (the third one is not shown for visualization 
purposes) are compared for the classical and effective solutions, it is clear that, far from the 
bounce, the classical and effective solutions agree. Near to the bounce we can see that the 
classical density diverges in a finite time while the effective one bounces. This also happens 
to the shear (that is not shown in the plot). The finiteness of these quantities is due to the 
bounce in the scale factors that now are not going to zero in a finite time (there are solutions 
where some scale factors do not bounce and continue approaching zero, but they need an 
infinite time to do so). This illustrates the manner in which the classical singularities are 
resolved. 

We understand by singularity resolution in the effective framework the possibility to 
evolve the solutions for an arbitrary time and that the solutions remain finite, thus signaling 
that the geodesies are inextendible. Classically, for homogeneous and anisotropic universes, 
the singularities are present when the scale factor goes to zero in a finite time. This is 
related with an infinite density (expansion and shear) and incompleteness of geodesies. If 
the density, expansion and shear remain finite and can be evolved for any time then we say 
that singularities are resolved, i.e., the scale factors are not equal zero in a finite time or, 
equivalently, the geodesies are inextendible [To] ITS]. 

The initial conditions for effective and classical solutions at t — 100 are: C\ = 0.01, C2 = 
0.02, c 3 = 0.03, pi = 10000,p 2 = 2000,p 3 = 200 (it is evolved back in time), the field 
momentum p^ = 108.9 is calculated from the Hamiltonian constraint. The initial conditions 
for the other classical solution at t — are: Ci = —1.920756 x 10 _3 ,c 2 = —7.971627 x 
lO- 2 , c 3 = -19.61028, pi = 15864.53, p 2 = 635.5782, p 3 = 4.317425 (and it is evolved forward 
in time). 

B. Isotropic Limit 

When x — > Bianchi II has Bianchi I as a limit. Furthermore, in this limit when 
Pi = P2 = P3 and ci = C2 = c 3 Bianchi I reduced to the isotropic case. Bianchi II solutions 
near the isotropic limit are shown in figure |2j where we can see that solutions to the effective 
equations have a maximal density equal to the critical density p cr i t [3 [T7j. The shear is 
close to zero (in this case it is less than 1.2 x 10~ 10 in Planck units) and presents a non 
trivial behavior because it has four maxima and vanishes when the density bounces. This is 
the first solution that is completely different from the known solutions in the isotropic and 
Bianchi I cases. The shear is small but not zero because this is not an isotropic solution, 
but it is only very close to it. One should keep in mind that the Bianchi I model and, 
therefore, the isotropic solutions are not contained within the Bianchi II solutions. 

In order to have control on the isotropy, we set the initial conditions at the bounce, i.e. 
when Ji\C\ = /12C2 = /23C3 = n/2 (if the solution is isotropic the three directions bounce at 
same time, but the opposite is not true) and then evolve back and forward in time, with 
these relations we only need three additional initial conditions, that can be either (ci, C2, C3) 
or (pi,P2;P3) or a combination of them. We take p\ — 1 x 10 5 , p<i = p 3 = 1 x 10 3 (which 
satisfy that x — > 0, i.e., p\ 3> P2P3) and using the fact that we are starting at the bounce 
(fiiCi = 7r/2) we calculate c\ = 2.185, C2 = 03 = 2.185 x 10 2 and from the Hamiltonian 
constraint p<f, = 2.86 x 10 5 . 
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FIG. 2: Isotropic limit. The density at the bounce is p/p C rit = 1; the shear is close to zero; the 
Kasner exponents are equal and change from 1/3 to 1/3 — 2/3 = —1/3 when they are evolved back 
in time; the Ricci scalar strongly changes its behavior near to the bounce and have its maximal 
value at the bounce. 

Note that in this case it is not true that pi = P2 = P3 and c\ = c 2 = C3. Then, one 
might why is this solution isotropic? The answer is that c\ and p\ are a rescaling of 02,3 and 
P2,3 that can be translated into a rescaling of the scalar factors such that the real criteria 
to say that it is isotropic are that the relations ai/a 2 , 02/03, a\ /a 3 remain constant and 
Hi = H 2 = H 3 , which are satisfied by our initial conditions. This can be imagined like an 
isotropic universe described by the evolution of a fiducial cuboid. The best way to see that 
this is an isotropic solution is to look at the Kasner exponents (Fig. [2]), which are the same 
ki = &2 = ^3 = ±1/3. The different signs specify when the directions are expanding (ki > 0) 
or contracting (ki < 0). 

C. Bianchi I Limit 

When x — > the Bianchi II model has Bianchi I as a limit, which has been extensively 
studied [10] and can be used as a reference point. In this limit we expect a maximal density 
less than the critical density p cr i t and a non zero shear with a maximal value at the bounce, 
also S 2 and Q must be conserved in the classical region (where S 2 + Q = 1, since K = 
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FIG. 3: Bianchi I limit. The density and shear have a dynamical contribution ($7 ^ 0, X 2 ^ 0); the 
shear reaches its maximal value (for bianchi I) at the bounce; the Kasner exponents change from 
ki to hi — 2/3 when they are evolved back in time, as was shown by |10| ; one direction is expanding 
(03), one contracting (02) and one bounces (ai); the expansion has one zero (bounce) and it is 
finite all the time; the Ricci scalar presents a slower changes in its behavior in comparison to the 
isotropic limit and it is finite too. 
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in Bianchi I). The Kasner exponents must satisfy the constraint equations for Bianchi I 
solutions (ki + k 2 + k 3 = ±1 and k\ + k 2 + k\ + k\ = 1) and must change like ki — > ki — 2/3 
when evolved back in time, as was shown in [TU]. All these facts are shown in figure [3j which 
tells us that in the Bianchi I limit, the effective dynamics of Bianchi II reproduces the know 
behavior for Bianchi I. 

A new and important feature of this solution is that it presents a maximal value of the 
shear for Bianchi I, as can be seen in figure [3j where the value of shear at the bounce is 
11.57 = 3° , as reported in [33j. Moreover, it can be noticed that just one direction a\ 
bounces and the other two directions a 2 , a 3 do not, i.e., H 2 , H 3 are not zero in a finite time. 
This implies that these directions continue going to zero (or infinite) but now they need an 
infinite time to reach these values. Note that in the classical region not all the directions 
are expanding (or contracting). This kind of universes have a classical singularity too, but 
it is a cigar-like singularity and it is different from the one in which all the directions are 
contracting (or expanding), called point-like singularity. Our simulations then show that this 
kind of singularity is resolved too, with the notion of 'singularity resolution' as explained in 
section IIV Al 

In this case the choice of initial conditions is: Ji\C\ = ir/2, fx 2 c 2 = tt/6, Ji 3 c 3 = 5ti/6 and 
Pi = 1 x 10 6 , p 2 = 100, p 3 = 100 (which imply that x = \/p2P3/pl ~ 0). With these initial 
conditions we get c x = 0.069, c 2 = 230.28, c 3 = 1151.39 and p^ = 5.84 x 10 4 . The initial 
time is at the global bounce (6 = 0), from which the solution is evolved back and forward 
in time. 



D. Locally Rotationally Symmetric (LRS) Solutions to Bianchi II 

These subclass of models are characterized by the directional scalar factors a±, a 2 , a 3 such 
that a 2 = a 3 . In particular, this means that at each point the space-time is invariant under 
rotation about a preferred spacelike axis. In our variables this is written as p 2 = p 3 and 
C2 = C3 (if these equalities are satisfied by the initial conditions then they are satisfied 
through the evolution). These solutions will be useful to study the limit when the density 
is equal to the maximal density. In Bianchi II one might expect a maximal density strictly 
less than 1.315p cr i t in all the solutions due to the presence of anisotropies. Here we show 
that this is not the case, and the density can achieve the maximal value as a consequence of 
the shear being zero at the bounce. This is very different from Bianchi I models where the 
presence of anisotropies makes the shear at the bounce always greater than zero. 

In order to have control on the density to make it maximal we put the initial conditions 
at the bounce, i.e., when \i\C\ = Ji 2 c 2 = Ji 3 c 3 = n/2 and then evolve back and forward in 
time. Also, we need to make x equal to the value that makes the density maximal, namely 
x = (i+7 2 )a ~ 0.83. If we take p 2 = p 3 = 1000 (LRS condition) then pi = 113 and using the 
fact that we are starting at the bounce (/2jC; = vr/2), we calculate c\ = 64.98, c 2 = c 3 = 7.34 
and, from the Hamiltonian constraint, p^ = 11029.97. 

The solution is plotted in figure |4j where it is shown that in fact the solution has a maximal 
density at the bounce. The shear is zero at the bounce and has a non trivial behavior, such 
as two maxima. In the evolution of the Kasner exponents we can see (looking from left 
to right) that two directions (a 2 ,a 3 ) are contracting and one direction (ai) is expanding, 
and the Kasner exponents look like Bianchi I exponents far away from the bounce. This is 
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FIG. 4: LRS solution. The density reaches its maximal value when the shear is zero at the bounce, 
but the shear is different from zero in the rest of the evolution; the dynamical observables at, Hi, kt 
for directions 2 and 3 are equal, this shows that it is a LRS solution; the directional scale factor 
a\ bounces one time and has two turn-around points; the curvature parameter K is different from 
zero and asymptotically it satisfies O + S 2 + K = 1 . 

because classically Bianchi II approaches Bianchi I when time goes to infinity [281 - 130] . then 
in the classical region the Kasner exponents of LRS Bianchi II have as a limit the Kasner 
exponents of LRS Bianchi I. In figure [4] we plot the Hubble parameter Hi, which is the 
one that presents a new behavior. From the plot it can be seen that the direction ai has 
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a bounce and two turn around points, which indicates also a new behavior that was not 
present in the isotropic and Bianchi I cases. The other two directions (a 2 , 03) bounces just 
one time for a total of three directional bounces and one global bounce (8 = 0). 

There is a valid question at this point. Why can the density reach its maximal value in 
Bianchi II and not in Bianchi I when there are anisotropies? This question arises because 
in both models there are new degrees of freedom due to the anisotropies, so in principle 
they will have a similar behavior respect to the distribution of the energy in gravitational 
waves. But the new feature in Bianchi II -not present in Bianchi I- is the nontrivial spatial 
curvature. This curvature gives a new degree of freedom with respect to the possible ways 
the energy density can be distributed. Now, the dynamical contribution not only comes 
from the matter density and the shear but also from the spatial curvature. This fact can 
be quantitatively understood from the fact that the curvature parameter K is non zero 
in Bianchi II. As we can see in figure [4j the plot for S 2 , K) shows that the curvature 
parameter K is different from zero. Thus, this provides also a qualitative explanation for 
the important difference between these two models, namely that the density can reach its 
maximal value at the bounce. This behavior must also occur in the Bianchi IX case [35] 
where there is spatial curvature as well. 



E. Vacuum Limit 

Finally we study the solutions in the vacuum limit with maximal shear. As in the previous 
limiting cases, we do not want to study the vacuum case (p = 0), but rather we want to 
study it as a limit approaching Bianchi II with a massless scalar field. That is, we shall 
study the solutions where p — > (or equivalently p^ — > 0). The problem is that the density 
goes to zero when the time goes to ±00 in all these solutions. Then we need to clarify what 
the vacuum limit means. In order to do this, we need to have in mind that the density has a 
maximal at the bounce, then we are interested in taking the vacuum limit p — > in a finite 
time (near to the bounce). That is, we want solutions with density near zero at the bounce 
or, equivalently, solutions where Q ~ asymptotically. Among these class of solutions we 
select those with maximal shear, this allows us to study the extreme solutions where all the 
dynamical contribution comes from the anisotropies. 

To obtain the initial conditions we found numerically the values of x and —n/2 < 
P2C2, < 37r/2 that makes the shear maximal with a density near to zero. To 
do this we use the analytical expression for the shear found in [31] and the Hamiltonian 



constraint Eq. (3.1). The values found are: x = 1.69, p\C\ = 0.5616, /2 2 c 2 = /2 3 c 3 = 1.94618. 
These values fix 4 initial conditions, namely pi,Ci,c 2 ,C3, as functions of the other two pa- 
rameters that are (in some sense) not relevant for the behavior that we want to study. In 
our simulations, we take p 2 = 100 and p% = 100. The initial time is at t = (vertical line in 
the plots). The solution is shown in figure [5j The first unexpected behavior is that the shear 
is zero at the bounce and reaches its maximal value on one side of the bounce. Moreover, p 2 
and p 3 have a new behavior because p 2 and p% are zero three times, these zeros can be seen 
in the plot for /ijCj, where /2 2 c 2 (and /2 3 c 3 ) cross the tt/2 line three times. This behavior is 
not present in Bianchi I nor in the solutions to Bianchi II studied in the previous sections. 
From the density plot we can see that it has a small value at the bounce and from the plot 
for the parameters (Q, T< 2 ,K), it can be appreciated that all the dynamical contribution 
comes from the anisotropies, as was expected. 
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FIG. 5: Vacuum limit. The shear reaches its maximal value at one side of the bounce and it is zero 
on the other side (the zero is not at the bounce) ; the density is not zero but has a small value (even 
at the bounce); piCi evolve from one classical solution (/ZjCj ~ 0) to other one (/2jQ pa 7r); p2 and p^ 
are zero three times; all the dynamical contributions come from the shear (fi « 0, K « 0,E 2 « 1); 
the Ricci scalar reaches its maximal value at one side of the bounce and a local minimum on the 
other side. 



19 



V. DISCUSSION 



In this paper, we analyzed the numerical solutions of the effective equations that come 
from the improved LQC dynamics of the Bianchi II model [12]. We choose a massless scalar 
field as the matter source. This effective theory comes from the construction of the full 
quantum theory and we expect that it gives some insights about the quantum dynamics 
of semiclassical states. The accuracy of the effective equations has been established in the 
isotropic cases and thus we expect that they should give an excellent approximation of the 
full quantum evolution for semiclassical states. 

Let us summarize our results. We considered the Bianchi II case at the classical and 
effective level. As Bianchi I is a limiting case for Bianchi II, we use the previous results for 
Bianchi I in order to have control on our model. As we have seen, we recover Bianchi I as 
a limiting case when x — > or equivalently p\ ^> P2P3- It is important to keep in mind that 
the Bianchi I model is a limiting case and is not contained within the Bianchi II model. 
The Bianchi I solutions are interesting by themselves, since they give information about the 
asymptotic behavior of Bianchi II. 

In order to determine how the classical singularities are resolved and how the effective 
equations evolve, we choose a set of observable quantities like density, shear, expansion, 
Ricci scalar, etc., and studied their evolution numerically. The equations of motion admit 
different limits that can be used to check the accuracy of the solutions and explore some new 
insights that the Bianchi II model offers. In order to systematically study these solutions, 
we started from the classical limit showing the way in which the effective solutions solve the 
singularities and reduce to the classical ones far away from the bounce. Next, we explored 
the isotropic limit included into the Bianchi I limit when there are no anisotropies. We 
found that these solutions have a maximal density equal to the critical density p cr it, the 
shear is close to zero and presents a non trivial behavior because it has four maxima and 
vanishes when the density bounces. The shear is not zero because this solution is not an 
exactly isotropic solution, but it is a solution very close to it. Recall that Bianchi I and, 
therefore, the isotropic solutions are not contained within the Bianchi II solutions. Later 
on, we added anisotropies to the Bianchi I limit and showed that they reproduce the known 
solutions [TD] . Here we have shown that the classical cigar-like singularities are resolved like 
the point-like singularities, with singularity resolution understood in terms of geometrical 
observables being well behaved. We could not show numerically the resolution of the barrel- 
like singularities because showing this implies a fine-tuning in the initial conditions, but we 
studied the limit of this kind of singularities, and there is nothing that indicates that they are 
not resolved, too. Then, we considered the Locally Rotationally Symmetric (LRS) model of 
Bianchi II and explored how to find the solutions with maximal density at the bounce. In this 
model at each point the space-time is invariant under rotation about a preferred spacelike 
axis. Here we showed that the density can have the maximal value as a consequence of the 
shear being zero at the bounce. This shear has a non trivial behavior, such as two maxima. 
The Kasner exponents look like Bianchi I exponents far away from the bounce, which is 
consistent with the fact that classically Bianchi II approaches Bianchi I when proper time 
goes to infinity [28-30J. It was then shown that one directional scale factor (ai) can change 
its behavior up to three times (it has one bounce and two turn-around points) and the other 
two directions (02, 03) bounce once, for a total of three directional bounces and one global 
bounce (when the expansion is zero). 

We also studied the solutions in the vacuum limit with maximal shear. This allowed 
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us to study the extreme solutions where all the dynamical contributions come from the 
anisotropies. These solutions present a small value of the density at the bounce and un- 
expected shear and Ricci scalar behaviors because they are asymmetrical, reaching their 
maximal value on one side of the bounce. We found that there are generic solutions in 
which the matter density is larger that its value in the isotropic solutions, with point-like 
and cigar-like singularities. Some important solutions are LRS, such as the one with maximal 
density and the vacuum limit. 

In order to have control over the solutions we found that the best way to do it is to put 
the initial conditions at the bounce (or near to it). There are also two important points that 
we have checked in this numerical work. The first one is the convergence of the solutions 
and the second one is the evolution of conserved quantities, additionally to the physical tests 
that the program must pass. 

Finally, these results can be used as a starting point to study the Bianchi IX model in 
order to know if the approach to the singularity of the effective solutions present Bianchi 
I behavior with Bianchi II transitions, as happens with the classical solutions in the BKL 
conjecture. In case of the full LQC dynamics it would be interesting to know whether the 
evolution of the semi-classical states reproduce all the new rich behavior that we get from 
the effective theory. From this point of view our work can be seen as the first step in this 
direction, since we already have a systematic study of the solutions of the effective theory. 
A similar study for the Bianchi IX model in underway |35j . 



We thank A. Ashtekar, Ed. Wilson-Ewing and D. Sloan for helpful discussions and com- 
ments. This work was in part supported by DGAPA-UNAM IN103610, by NSF PHY0854743 
and by the Eberly Research Funds of Penn State. 



There are two important points that need to be cover when there is a numerical work: 
one is the convergence of solutions and other is the evolution of conserved quantities. This 
means that, additionally to the physical tests that the program must pass, also must present 
a convergence of solutions, i.e., numerical solutions approach to analytical ones when the 
accuracy is improved, this say us that we are near to the analytical solution with a small 
relative error. Remember that numerical solutions (in general) are never on the analytical 
solutions, all we can say is that they converge to them. Numerical solutions must also evolve 
on the constraint surface and they must preserve the conserved quantities, this ensures that 
they are evolving on the physical phase space. In figure [6] it is shown the convergence 
of Hamiltonian constraint (Ch ~ 0) when the time step is reduced, this also show that 
the constraint (or equivalently p^) is conserved. The quantities plotted in figure [6] are the 
relative error for the constraint 
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Appendix A: Convergence and Conservation 
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with different resolutions, and the error functions 
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where the subindices in the Hamiltonian constraint mean resolution 2 (with d^) and 
resolution 1 (with dti), where dt 2 = dt 1 /2. The method used to integrate the equa- 
tions is a Runge-Kutta 4 (RK4), while the resolutions used for the convergence tests are 
dt = 0.01,0.005,0.0025,0.00125,0.000625. The error functions for (c u pi) present similar 
behaviors. We can define the convergence order as 



n 



h-f2 

J 2 — f3 



(A3) 



where fi is any evolved function at resolution i, with dti > dti+i. The convergence factor n 
for a RK4 must be n = 2 4 = 16. We obtain in our solutions n w 16.2, which say us that 
solutions convergence as fast as expected. 



Convergence of the Constraint 




FIG. 6: Convergence of the constraint. The first plot shows the convergence of the solutions and 
that the constraint is conserved. The second plot shows the error functions versus the bigger 
resultion used to calculate their. 
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